# Colección de paquetes de Tidyverse
library(tidyverse)
# Estilos para ggplot2
library(ggthemes)
# Paletas de colores de RColorBrewer
library(RColorBrewer)
# Paletas de colores de viridis
library(viridisLite)
# Gráficos interactivos
library(plotly)
# Manejo de datos vectoriales
library(sf)
# Manejo de datos raster
library(terra)
# Manejo de datos raster
library(raster)
# Mapas interactivos
library(leaflet)
# Acceso a datos en GBIF
library(rgbif)
# Datos geoespaciales
library(geodata)
# Modelado de distribución de especies
library(dismo)
library(rJava)Modelos de Nichos Ecológicos
Carga de paquetes
Obtención de datos de presencia
# Nombre de la especie
especie <- "Bradypus variegatus"# Consulta a GBIF
respuesta <- occ_search(
scientificName = especie,
hasCoordinate = TRUE,
hasGeospatialIssue = FALSE,
limit = 10000
)
# Extraer datos de presencia
presencia <- respuesta$data# Guardar los datos de presencia en un archivo CSV
write_csv(presencia, 'presencia.csv')# Leer los datos de presencia de un archivo CSV
presencia <- read_csv('presencia.csv')Warning: One or more parsing issues, call `problems()` on your data frame for details,
e.g.:
dat <- vroom(...)
problems(dat)
Rows: 6619 Columns: 173
── Column specification ────────────────────────────────────────────────────────
Delimiter: ","
chr (122): scientificName, issues, datasetKey, publishingOrgKey, installati...
dbl (28): key, decimalLatitude, decimalLongitude, crawlId, taxonKey, kingd...
lgl (17): isSequenced, isInCluster, acceptedNameUsage, locationRemarks, ta...
dttm (5): lastCrawled, lastParsed, dateIdentified, modified, lastInterpreted
date (1): georeferencedDate
ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
# Convertir a conjunto de datos espacial
presencia <- st_as_sf(
presencia,
coords = c("decimalLongitude", "decimalLatitude"),
remove = FALSE, # conservar las columnas de las coordenadas
crs = 4326
)# Gráfico ggplot2
grafico_ggplot2 <-
presencia |>
st_drop_geometry() |>
group_by(year) |>
summarize(n = n()) |>
ggplot(aes(x = year, y = n)) +
geom_line() +
geom_point(
aes(
text = paste0(
"Año: ", year, "\n",
"Cantidad de registros: ", n
)
)
) +
ggtitle("Cantidad de registros de presencia por año") +
xlab("Año") +
ylab("Cantidad de registros de presencia") +
labs(caption = "Fuente: GBIF") +
theme_economist()Warning in geom_point(aes(text = paste0("Año: ", year, "\n", "Cantidad de
registros: ", : Ignoring unknown aesthetics: text
# Gráfico plotly
ggplotly(grafico_ggplot2, tooltip = "text") |>
config(locale = 'es')# Gráfico ggplot2
grafico_ggplot2 <-
presencia |>
st_drop_geometry() |>
group_by(year) |>
summarize(n = n()) |>
ggplot(aes(x = year, y = n)) +
geom_line() +
geom_point(
aes(
text = paste0(
"Año: ", year, "\n",
"Cantidad de registros: ", n
)
)
) +
ggtitle("Cantidad de registros de presencia por año") +
xlab("Año") +
ylab("Cantidad de registros de presencia") +
labs(caption = "Fuente: GBIF") +
theme_economist()Warning in geom_point(aes(text = paste0("Año: ", year, "\n", "Cantidad de
registros: ", : Ignoring unknown aesthetics: text
# Gráfico plotly
ggplotly(grafico_ggplot2, tooltip = "text") |>
config(locale = 'es')# Mapa
leaflet() |>
addTiles(group = "Mapa general") |>
addProviderTiles(
providers$Esri.WorldImagery,
group = "Imágenes satelitales"
) |>
addProviderTiles(
providers$CartoDB.Positron,
group = "Mapa blanco"
) |>
addCircleMarkers(
# capa de registros de presencia (puntos)
data = presencia,
stroke = F,
radius = 3,
fillColor = 'red',
fillOpacity = 1,
popup = paste(
paste0("<strong>País: </strong>", presencia$country),
paste0("<strong>Localidad: </strong>", presencia$locality),
paste0("<strong>Fecha: </strong>", presencia$eventDate),
paste0("<strong>Fuente: </strong>", presencia$institutionCode),
paste0("<a href='", presencia$occurrenceID, "'>Más información</a>"),
sep = '<br/>'
),
group = "Registros de Bradypus variegatus"
) |>
addLayersControl(
baseGroups = c("Mapa general", "Imágenes satelitales", "Mapa blanco"),
overlayGroups = c("Registros de Bradypus variegatus"))# Consulta a WorldClim
clima <- worldclim_global(var = 'bio', res = 10, path = tempdir())
# Nombres de las variables climáticas
names(clima) [1] "wc2.1_10m_bio_1" "wc2.1_10m_bio_2" "wc2.1_10m_bio_3" "wc2.1_10m_bio_4"
[5] "wc2.1_10m_bio_5" "wc2.1_10m_bio_6" "wc2.1_10m_bio_7" "wc2.1_10m_bio_8"
[9] "wc2.1_10m_bio_9" "wc2.1_10m_bio_10" "wc2.1_10m_bio_11" "wc2.1_10m_bio_12"
[13] "wc2.1_10m_bio_13" "wc2.1_10m_bio_14" "wc2.1_10m_bio_15" "wc2.1_10m_bio_16"
[17] "wc2.1_10m_bio_17" "wc2.1_10m_bio_18" "wc2.1_10m_bio_19"
# Definir la extensión del área de estudio
area_estudio <- ext(
min(presencia$decimalLongitude) - 5,
max(presencia$decimalLongitude) + 5,
min(presencia$decimalLatitude) - 5,
max(presencia$decimalLatitude) + 5
)
# Recortar las variables bioclimáticas al área de estudio
clima <- crop(clima, area_estudio)# Paleta de colores de temperatura
colores_temperatura <- colorNumeric(
# palette = "inferno",
# palette = "magma",
palette = rev(brewer.pal(11, "RdYlBu")),
values(clima$wc2.1_10m_bio_1),
na.color = "transparent"
)
# Paleta de colores de precipitación
colores_precipitacion <- colorNumeric(
# palette = "viridis",
# palette = "YlGnBu",
palette = "Blues",
values(clima$wc2.1_10m_bio_12),
na.color = "transparent"
)
# Mapa
leaflet() |>
addTiles(group = "Mapa general") |>
addProviderTiles(
providers$Esri.WorldImagery,
group = "Imágenes satelitales"
) |>
addProviderTiles(
providers$CartoDB.Positron,
group = "Mapa blanco"
) |>
addRasterImage( # capa raster de temperatura
clima$wc2.1_10m_bio_1,
colors = colores_temperatura, # paleta de colores
opacity = 0.6,
group = "Temperatura",
) |>
addRasterImage( # capa raster de precipitación
clima$wc2.1_10m_bio_12,
colors = colores_precipitacion, # paleta de colores
opacity = 0.6,
group = "Precipitación",
) |>
addCircleMarkers(
# capa de registros de presencia (puntos)
data = presencia,
stroke = F,
radius = 3,
fillColor = 'red',
fillOpacity = 1,
popup = paste(
paste0("<strong>País: </strong>", presencia$country),
paste0("<strong>Localidad: </strong>", presencia$locality),
paste0("<strong>Fecha: </strong>", presencia$eventDate),
paste0("<strong>Fuente: </strong>", presencia$institutionCode),
paste0("<a href='", presencia$occurrenceID, "'>Más información</a>"),
sep = '<br/>'
),
group = "Registros de Bradypus variegatus"
) |>
addLegend(
title = "Temperatura",
values = values(clima$wc2.1_10m_bio_1),
pal = colores_temperatura,
position = "bottomleft",
group = "Temperatura"
) |>
addLegend(
title = "Precipitación",
values = values(clima$wc2.1_10m_bio_12),
pal = colores_precipitacion,
position = "bottomleft",
group = "Precipitación"
) |>
addLayersControl(
# control de capas
baseGroups = c("Mapa general", "Imágenes satelitales", "Mapa blanco"),
overlayGroups = c("Temperatura", "Precipitación", "Registros de Bradypus variegatus")
) |>
hideGroup("Precipitación")coordenadas_presencia <- data.frame(
decimalLongitude = presencia$decimalLongitude,
decimalLatitude = presencia$decimalLatitude
)
coordenadas_presencia <- unique(coordenadas_presencia)# Dividir los datos en entrenamiento y evaluación
set.seed(123)
n_presencia <- nrow(coordenadas_presencia)
indices_entrenamiento <- sample(
1:n_presencia,
size = round(0.7 * n_presencia)
)
entrenamiento <- coordenadas_presencia[indices_entrenamiento, ]
evaluacion <- coordenadas_presencia[-indices_entrenamiento, ]# Los datos de clima deben convertirse al formato que usa el paquete raster
# debido a es este el que acepta el paquete dismo
clima_raster <- raster::stack(clima)
# Ejecutar el modelo
#modelo_maxent <- maxent(x = clima_raster, p = entrenamiento)
modelo_maxent <- domain(x = clima_raster, p = entrenamiento)
# Crear un raster basado en el modelo
prediccion <- predict(modelo_maxent, clima_raster)# Extraer valores del raster de predicción para los puntos de evaluación de presencias
eval_pres <- terra::extract(
prediccion,
evaluacion[, c('decimalLongitude', 'decimalLatitude')]
)
# Generar puntos aleatorios de ausencia
ausencias <- randomPoints(mask = clima_raster, n = 1000)
# Extraer valores del raster de predicción para los puntos de ausencias
eval_aus <- terra::extract(
prediccion,
ausencias
)
# Generar estadísticas de evaluación del modelo
evaluation <- evaluate(p = eval_pres, a = eval_aus)
# Graficación de la curva ROC
plot(evaluation, 'ROC')
# Paleta de colores del modelo
colores_modelo <- colorNumeric(
palette = c("white", "black"),
values(prediccion),
na.color = "transparent"
)
# Mapa
leaflet() |>
addTiles(group = "Mapa general") |>
addProviderTiles(
providers$Esri.WorldImagery,
group = "Imágenes satelitales"
) |>
addProviderTiles(
providers$CartoDB.Positron,
group = "Mapa blanco"
) |>
addRasterImage( # capa raster de temperatura
clima$wc2.1_10m_bio_1,
colors = colores_temperatura, # paleta de colores
opacity = 0.6,
group = "Temperatura",
) |>
addRasterImage( # capa raster de precipitación
clima$wc2.1_10m_bio_12,
colors = colores_precipitacion, # paleta de colores
opacity = 0.6,
group = "Precipitación",
) |>
addRasterImage( # capa raster del modelo de distribución
prediccion,
colors = colores_modelo,
opacity = 0.6,
group = "Modelo de distribución",
) |>
addCircleMarkers(
# capa de registros de presencia (puntos)
data = presencia,
stroke = F,
radius = 3,
fillColor = 'red',
fillOpacity = 1,
popup = paste(
paste0("<strong>País: </strong>", presencia$country),
paste0("<strong>Localidad: </strong>", presencia$locality),
paste0("<strong>Fecha: </strong>", presencia$eventDate),
paste0("<strong>Fuente: </strong>", presencia$institutionCode),
paste0("<a href='", presencia$occurrenceID, "'>Más información</a>"),
sep = '<br/>'
),
group = "Registros de Bradypus variegatus"
) |>
addLegend(
title = "Temperatura",
values = values(clima$wc2.1_10m_bio_1),
pal = colores_temperatura,
position = "bottomleft",
group = "Temperatura"
) |>
addLegend(
title = "Precipitación",
values = values(clima$wc2.1_10m_bio_12),
pal = colores_precipitacion,
position = "bottomleft",
group = "Precipitación"
) |>
addLegend(
title = "Modelo de distribución",
values = values(prediccion),
pal = colores_modelo,
position = "bottomright",
group = "Modelo de distribución"
) |>
addLayersControl(
# control de capas
baseGroups = c("Mapa general", "Imágenes satelitales", "Mapa blanco"),
overlayGroups = c(
"Temperatura",
"Precipitación",
"Modelo de distribución",
"Registros de Bradypus variegatus"
)
) |>
hideGroup("Temperatura") |>
hideGroup("Precipitación")Warning in colors(.): Some values were outside the color scale and will be
treated as NA